      subroutine initial_state(np,kr,harris,t_shape,flvx,flvy,flvz,bfldx,bfldy,bfldz,efldx,efldy,efldz)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      use vast_kind_param, ONLY:  double
      use corgan_com_M
      use cindex_com_M
      use numpar_com_M
      use cprplt_com_M
      use cophys_com_M
      use ctemp_com_M
      use blcom_com_M
!
      implicit none
!-----------------------------------------------
!   G l o b a l   P a r a m e t e r s
!-----------------------------------------------
     integer :: np,kr
     real(double) :: harris,flvx,flvy,flvz,bfldx,bfldy,bfldz,t_shape,efldx,efldy,efldz
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
     real(double) :: xc,zc,zone,ztwo
!-----------------------------------------------
!   E x t e r n a l   F u n c t i o n s
!!-----------------------------------------------
!-----------------------------------------------
!
!
!     CASES
!       1- Harris (Tail) generalized by Birn for tail
!       2- Force Free by Lapenta&Delzanno
!       3- Fadeev by Wan
!       4- Sheared flow and magnetic field
!       5- Rope equilibrium
!       6- Equilibrium loaded from file + uniform
!       7- Equilibrium uniform
!       8- One wave equilibrium
!	    9- expanding bubble
!
! density perturbation that is consistent with initial insland perturbation
!                        harris = harris - pert_gem*bx_ini*((2.d0*pi/(xr-xl))**2 + &
!                             (pi/(zf-ze))**2)*cos(2.d0*pi*(xc-0.5d0*(xr+xl))/(xr-xl)) &
!                             *cos(pi*(zc-.5d0*(zf+ze))/(zf-ze))
flvx = 1.
flvy = 1.
flvz = 1.
efldx=0.
efldy=0.
efldz=0.
!write(*,*)'initial_state',initial_equilibrium
select case(initial_equilibrium)
    case(1)
        call tail_equilibrium(np,kr,harris,t_shape,flvx,flvy,flvz,bfldx,bfldy,bfldz)
    case(2)
        call forcefree_equilibrium(np,kr,harris,t_shape,flvx,flvy,flvz,bfldx,bfldy,bfldz)
    case(3)
        call fadeev_equilibrium(np,kr,harris,t_shape,flvx,flvy,flvz,bfldx,bfldy,bfldz)
    case(4)
        call shear_equilibrium(np,kr,harris,t_shape,flvx,flvy,flvz,bfldx,bfldy,bfldz)
    case(5)
        call rope_equilibrium(np,kr,harris,t_shape,flvx,flvy,flvz,bfldx,bfldy,bfldz)
    case(6)
        call uniform_equilibrium(np,kr,harris,t_shape,flvx,flvy,flvz,bfldx,bfldy,bfldz)
    case(7)
        call uniform_equilibrium(np,kr,harris,t_shape,flvx,flvy,flvz,bfldx,bfldy,bfldz)
    case(8)
        call onewave_equilibrium(np,kr,harris,t_shape,flvx,flvy,flvz,bfldx,bfldy,bfldz,efldx,efldy,efldz)
    case(9)
        call bubble_equilibrium(np,kr,harris,t_shape,flvx,flvy,flvz,bfldx,bfldy,bfldz)
    case(10)
        call doubletearing_equilibrium(np,kr,harris,t_shape,flvx,flvy,flvz,bfldx,bfldy,bfldz)
    case default
        write(*,*)'The initial state is undefined, check the input file and try again'
        stop
end select
!
end subroutine initial_state

subroutine tail_equilibrium(np,kr,harris,t_shape,flvx,flvy,flvz,bfldx,bfldy,bfldz)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param
      use corgan_com_M
      use cindex_com_M
      use numpar_com_M
      use cprplt_com_M
      use cophys_com_M
      use ctemp_com_M
      use blcom_com_M
!
      implicit none
!-----------------------------------------------
!   G l o b a l   P a r a m e t e r s
!-----------------------------------------------
     real(double) :: harris,flvx,flvy,flvz,bfldx,bfldy,bfldz,t_shape
     integer :: np,kr
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
     real(double)::fx,fxone,fxtwo,segno,xc,zc,zone,ztwo,xone,xtwo
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
!-----------------------------------------------
!   E x t e r n a l   F u n c t i o n s
!-----------------------------------------------
!-----------------------------------------------!
!    Compute grid quantities
!
!
                     ipjk = ijk + iwid
                     ipjpk = ijk + iwid + jwid
                     ijpk = ijk + jwid
!
                     ijkp = ijk + kwid
                     ipjkp = ijk + iwid + kwid
                     ijpkp = ijk + jwid + kwid
                     ipjpkp = ijk + iwid + jwid + kwid
!

                     xc = 0.125d0*(x(ipjk)+x(ipjpk)+x(ijpk)+x(ijk)+x(ipjkp)+x(ipjpkp)+x&
                          (ijpkp)+x(ijkp))
                     zc = 0.125d0*(z(ipjk)+z(ipjpk)+z(ijpk)+z(ijk)+z(ipjkp)+z(ipjpkp)+z&
                          (ijpkp)+z(ijkp))
!
                     zone = 0.25d0*(z(ipjk)+z(ipjpk)+z(ijpk)+z(ijk))
                     ztwo = 0.25d0*(z(ipjkp)+z(ipjpkp)+z(ijpkp)+z(ijkp))
                     xtwo = 0.25d0*(x(ipjk)+x(ipjpk)+x(ipjkp)+x(ipjpkp))
                     xone = 0.25d0*(x(ijk)+x(ijpk)+x(ijkp)+x(ijpkp))
!


!       BIRN profile
!
fx = (1. + eps_rec/rnu_rec*px(np)/el_rec)**(-rnu_rec)
fxone = (1. + eps_rec/rnu_rec*xone/el_rec)**(-rnu_rec)
fxtwo = (1. + eps_rec/rnu_rec*xtwo/el_rec)**(-rnu_rec)
!
!       Quasiparabolic profile (Pritchett)
!
!        fx=exp(-px(np)*eps_rec/el_rec)
!        fxone=exp(-xone*eps_rec/el_rec)
!        fxtwo=exp(-xtwo*eps_rec/el_rec)

!
segno = sign(1.,qom(icoi(kr)))

if (abs(vvi(kr)) < 1.E-10) then
   harris = 1.
else
   harris = el_rec*fx/(ztwo - zone)*(tanh((ztwo - 0.5*(zf&
            + ze))*fx/el_rec) - tanh((zone - 0.5*(zf + ze))*fx/el_rec))
end if

t_shape=1d0

flvx=uvi(kr)*flvx*segno + uvip(kr)*flvx*segno*cos(2.*pi/(yt - yb)*mpert*py(np))
flvy=vvi(kr)*flvy*segno + vvip(kr)*flvx*segno*cos(2.*pi/(yt - yb)*mpert*py(np))
flvz=wvi(kr)*flvz*segno + wvip(kr)*flvx*segno*cos(2.*pi/(yt - yb)*mpert*py(np))

! campo magnetico
!     bfldx=bx_ini*tanh((zc-.5*(zf+ze))/el_rec)
bfldx = bx_ini*el_rec*(log(cosh((ztwo - 0.5*(zf + ze))*fx/&
        el_rec)/fx) - log(cosh((zone - 0.5*(zf + ze))*fx/el_rec)/fx))&
        /(ztwo - zone)
bfldy = by_ini
!bfldy = by_ini*el_rec*(log(cosh((ztwo - 0.5*(zf + ze))*fx/&
!        el_rec)/fx) - log(cosh((zone - 0.5*(zf + ze))*fx/el_rec)/fx))&
!        /(ztwo - zone)
!               bfldy = by_ini*el_rec*(log(cosh((ztwo - 0.5*(zf + ze))/el_rec&
!                  )) - log(cosh((zone - 0.5*(zf + ze))/el_rec)))/(ztwo - zone)
bfldz = (-bz_ini) - bx_ini*el_rec*(log(cosh((zc - 0.5*(zf + &
        ze))*fxtwo/el_rec)/fxtwo) - log(cosh((zc - 0.5*(zf + ze))*&
        fxone/el_rec)/fxone))/(xtwo - xone)

end subroutine tail_equilibrium

subroutine forcefree_equilibrium(np,kr,harris,t_shape,flvx,flvy,flvz,bfldx,bfldy,bfldz)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param
      use corgan_com_M
      use cindex_com_M
      use numpar_com_M
      use cprplt_com_M
      use cophys_com_M
      use ctemp_com_M
      use blcom_com_M
!
      implicit none
!-----------------------------------------------
!   G l o b a l   P a r a m e t e r s
!-----------------------------------------------
     real(double) :: harris,flvx,flvy,flvz,bfldx,bfldy,bfldz,t_shape
     integer :: np,kr
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
   real(double) :: fx,fxone,fxtwo,shaperx,shapery,segno,xc,zc,xone,xtwo,zone,ztwo
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
!-----------------------------------------------
!   E x t e r n a l   F u n c t i o n s
!-----------------------------------------------
!-----------------------------------------------!
!    Compute grid quantities
!
!
                     ipjk = ijk + iwid
                     ipjpk = ijk + iwid + jwid
                     ijpk = ijk + jwid
!
                     ijkp = ijk + kwid
                     ipjkp = ijk + iwid + kwid
                     ijpkp = ijk + jwid + kwid
                     ipjpkp = ijk + iwid + jwid + kwid
!

                     xc = 0.125d0*(x(ipjk)+x(ipjpk)+x(ijpk)+x(ijk)+x(ipjkp)+x(ipjpkp)+x&
                          (ijpkp)+x(ijkp))
                     zc = 0.125d0*(z(ipjk)+z(ipjpk)+z(ijpk)+z(ijk)+z(ipjkp)+z(ipjpkp)+z&
                          (ijpkp)+z(ijkp))
!
                     zone = 0.25d0*(z(ipjk)+z(ipjpk)+z(ijpk)+z(ijk))
                     ztwo = 0.25d0*(z(ipjkp)+z(ipjpkp)+z(ijpkp)+z(ijkp))
                     xtwo = 0.25d0*(x(ipjk)+x(ipjpk)+x(ipjkp)+x(ipjpkp))
                     xone = 0.25d0*(x(ijk)+x(ijpk)+x(ijkp)+x(ijpkp))
!
harris=1d0
t_shape=1d0
fx=1.d0
fxone=1d0
fxtwo=1d0
segno = sign(1.,qom(icoi(kr)))
shaperx = tanh((pz(np)-0.5*(zf+ze))/el_rec)/cosh((pz(np)-0.5*(zf+ze))/el_rec)/el_rec
shapery = 1.d0/cosh((pz(np)-0.5*(zf+ze))/el_rec)**2/el_rec
!
flvx=uvi(kr)*flvx*segno*shaperx+ uvip(kr)*flvx*segno*cos(2.*pi/(yt - yb)*mpert*py(np))
flvy=vvi(kr)*flvy*segno*shapery+ vvip(kr)*flvx*segno*cos(2.*pi/(yt - yb)*mpert*py(np))
flvz=wvi(kr)*flvz*segno+ wvip(kr)*flvx*segno*cos(2.*pi/(yt - yb)*mpert*py(np))
! uvi,vvi must be chosen to support the current
! campo magnetico

bfldx = bx_ini*el_rec*(log(cosh((ztwo - 0.5*(zf + ze))*fx/&
        el_rec)/fx) - log(cosh((zone - 0.5*(zf + ze))*fx/el_rec)/fx))&
        /(ztwo - zone)
bfldy = by_ini/cosh((zc-0.5*(zf + ze))/el_rec)
bfldz = bz_ini

end subroutine forcefree_equilibrium

subroutine fadeev_equilibrium(np,kr,harris,t_shape,flvx,flvy,flvz,bfldx,bfldy,bfldz)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param
      use corgan_com_M
      use cindex_com_M
      use numpar_com_M
      use cprplt_com_M
      use cophys_com_M
      use ctemp_com_M
      use blcom_com_M
!
      implicit none
!-----------------------------------------------
!   G l o b a l   P a r a m e t e r s
!-----------------------------------------------
     real(double) :: harris,flvx,flvy,flvz,bfldx,bfldy,bfldz,t_shape
     integer :: np,kr
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
real(double) :: segno,xc,zc
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
!-----------------------------------------------
!   E x t e r n a l   F u n c t i o n s
!-----------------------------------------------
!-----------------------------------------------!
!    Compute grid quantities
!
!
                     ipjk = ijk + iwid
                     ipjpk = ijk + iwid + jwid
                     ijpk = ijk + jwid
!
                     ijkp = ijk + kwid
                     ipjkp = ijk + iwid + kwid
                     ijpkp = ijk + jwid + kwid
                     ipjpkp = ijk + iwid + jwid + kwid
!

                     xc = 0.125d0*(x(ipjk)+x(ipjpk)+x(ijpk)+x(ijk)+x(ipjkp)+x(ipjpkp)+x&
                          (ijpkp)+x(ijkp))
                     zc = 0.125d0*(z(ipjk)+z(ipjpk)+z(ijpk)+z(ijk)+z(ipjkp)+z(ipjpkp)+z&
                          (ijpkp)+z(ijkp))

!    Fadeev island chain equilibrium profile
segno = sign(1.,qom(icoi(kr)))
harris = (1.d0-epsil**2)/ (epsil*cos((xc-.5*(xr+xl))/el_rec) + cosh((zc-.5*(zf+ze))/el_rec))**2

flvx=uvi(kr)*flvx*segno + uvip(kr)*flvx*segno*cos(2.*pi/(yt - yb)*mpert*py(np))
flvy=vvi(kr)*flvy*segno + vvip(kr)*flvx*segno*cos(2.*pi/(yt - yb)*mpert*py(np))
flvz=wvi(kr)*flvz*segno + wvip(kr)*flvx*segno*cos(2.*pi/(yt - yb)*mpert*py(np))

t_shape=1d0
! campo magnetico

bfldx = bx_ini*sinh((zc-.5*(zf+ze))/el_rec)/ &
        (epsil*cos((xc-.5*(xr+xl))/el_rec) + cosh((zc-.5*(zf+ze))/el_rec))
bfldy = by_ini
bfldz = epsil*bx_ini*sin((xc-.5*(xr+xl))/el_rec)/ &
        (epsil*cos((xc-.5*(xr+xl))/el_rec) + cosh((zc-.5*(zf+ze))/el_rec))

end subroutine fadeev_equilibrium

subroutine shear_equilibrium(np,kr,harris,t_shape,flvx,flvy,flvz,bfldx,bfldy,bfldz)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param
      use corgan_com_M
      use cindex_com_M
      use numpar_com_M
      use cprplt_com_M
      use cophys_com_M
      use ctemp_com_M
      use blcom_com_M
!
      implicit none
!-----------------------------------------------
!   G l o b a l   P a r a m e t e r s
!-----------------------------------------------
     real(double) :: harris,flvx,flvy,flvz,bfldx,bfldy,bfldz,t_shape
     integer :: np,kr
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
     real(double)::fx,shaper,segno,xc,zc,zone,ztwo
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
!-----------------------------------------------
!   E x t e r n a l   F u n c t i o n s
!-----------------------------------------------
!
!    Compute grid quantities
!
!
                     ipjk = ijk + iwid
                     ipjpk = ijk + iwid + jwid
                     ijpk = ijk + jwid
!
                     ijkp = ijk + kwid
                     ipjkp = ijk + iwid + kwid
                     ijpkp = ijk + jwid + kwid
                     ipjpkp = ijk + iwid + jwid + kwid
!

                     xc = 0.125d0*(x(ipjk)+x(ipjpk)+x(ijpk)+x(ijk)+x(ipjkp)+x(ipjpkp)+x&
                          (ijpkp)+x(ijkp))
                     zc = 0.125d0*(z(ipjk)+z(ipjpk)+z(ijpk)+z(ijk)+z(ipjkp)+z(ipjpkp)+z&
                          (ijpkp)+z(ijkp))
!
                     zone = 0.25d0*(z(ipjk)+z(ipjpk)+z(ijpk)+z(ijk))
                     ztwo = 0.25d0*(z(ipjkp)+z(ipjpkp)+z(ijpkp)+z(ijkp))
!
!       BIRN profile
!
fx=1d0

segno = sign(1.,qom(icoi(kr)))
harris = el_rec*fx/(ztwo - zone)*(tanh((ztwo - 0.5*(zf&
              + ze))*fx/el_rec) - tanh((zone - 0.5*(zf + ze))*fx/el_rec))

t_shape=1d0

shaper = tanh((pz(np)-0.5*(zf+ze))/el_rec)

flvx=uvi(kr)*flvx*segno + uvip(kr)/3.*flvx*segno*cos(2.*pi/(xr - xl)*mpert*px(np)) + uvip(kr)*shaper
flvy=vvi(kr)*flvy*segno + vvip(kr)/3.*flvx*segno*cos(2.*pi/(xr - xl)*mpert*px(np)) + vvip(kr)*shaper
flvz=wvi(kr)*flvz*segno + uvip(kr)/3.*flvx*segno*sin(2.*pi/(xr - xl)*mpert*px(np)) + wvip(kr)*shaper
! campo magnetico

bfldx = bx_ini*el_rec*(log(cosh((ztwo - 0.5*(zf + ze))*fx/&
        el_rec)/fx) - log(cosh((zone - 0.5*(zf + ze))*fx/el_rec)/fx))&
        /(ztwo - zone)
bfldy = by_ini
bfldz =-bz_ini


end subroutine shear_equilibrium

subroutine uniform_equilibrium(np,kr,harris,t_shape,flvx,flvy,flvz,bfldx,bfldy,bfldz)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param
      use corgan_com_M
      use cindex_com_M
      use numpar_com_M
      use cprplt_com_M
      use cophys_com_M
      use ctemp_com_M
      use blcom_com_M
!
      implicit none
!-----------------------------------------------
!   G l o b a l   P a r a m e t e r s
!-----------------------------------------------
     real(double) :: harris,flvx,flvy,flvz,bfldx,bfldy,bfldz,t_shape
     integer :: np,kr
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
!-----------------------------------------------
!   E x t e r n a l   F u n c t i o n s
!-----------------------------------------------
!-----------------------------------------------!
! Unifrom density
harris=1d0
! Uniform temperature
t_shape=1.d0
!
! Plasma flow
!
flvx=uvi(kr)!+ uvip(kr)*cos(2.*pi/(kbar*dz)*mpert*pz(np))
flvy=vvi(kr)!+ vvip(kr)*cos(2.*pi/(kbar*dz)*mpert*pz(np))
flvz=wvi(kr)!+ wvip(kr)*cos(2.*pi/(kbar*dz)*mpert*pz(np))
!
! Magnetic Field
!     bfldx=bx_ini*tanh((zc-.5*(zf+ze))/el_rec)
bfldx = bx_ini
bfldy = by_ini
bfldz = bz_ini

end subroutine uniform_equilibrium

subroutine bubble_equilibrium(np,kr,harris,t_shape,flvx,flvy,flvz,bfldx,bfldy,bfldz)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param
      use corgan_com_M
      use cindex_com_M
      use numpar_com_M
      use cprplt_com_M
      use cophys_com_M
      use ctemp_com_M
      use blcom_com_M
!
      implicit none
!-----------------------------------------------
!   G l o b a l   P a r a m e t e r s
!-----------------------------------------------
     real(double) :: harris,flvx,flvy,flvz,bfldx,bfldy,bfldz,t_shape
     integer :: np,kr
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
     real(double)::xc,yc,zc,rc
!-----------------------------------------------
!   E x t e r n a l   F u n c t i o n s
!-----------------------------------------------
!-----------------------------------------------!
!
!    Compute grid quantities
!
!
                     ipjk = ijk + iwid
                     ipjpk = ijk + iwid + jwid
                     ijpk = ijk + jwid
!
                     ijkp = ijk + kwid
                     ipjkp = ijk + iwid + kwid
                     ijpkp = ijk + jwid + kwid
                     ipjpkp = ijk + iwid + jwid + kwid
!

                     xc = 0.125d0*(x(ipjk)+x(ipjpk)+x(ijpk)+x(ijk)+x(ipjkp)+x(ipjpkp)+x&
                          (ijpkp)+x(ijkp))
                     yc = 0.125d0*(y(ipjk)+y(ipjpk)+y(ijpk)+y(ijk)+y(ipjkp)+y(ipjpkp)+y&
                          (ijpkp)+y(ijkp))
                     zc = 0.125d0*(z(ipjk)+z(ipjpk)+z(ijpk)+z(ijk)+z(ipjkp)+z(ipjpkp)+z&
                          (ijpkp)+z(ijkp))
                     rc = sqrt(xc**2+yc**2+zc**2)
!
! density
if(rc<el_rec) then
harris=1d0
t_shape=1d0
else
harris=.1d0
t_shape=.1d0
end if

!
! Plasma flow
!
flvx=uvi(kr)
flvy=vvi(kr)
flvz=wvi(kr)
!
! Magnetic Field
!     bfldx=bx_ini*tanh((zc-.5*(zf+ze))/el_rec)
bfldx = bx_ini
bfldy = by_ini
bfldz = bz_ini

end subroutine bubble_equilibrium

subroutine rope_equilibrium(np,kr,harris,vt_shape,flvx,flvy,flvz,bfldx,bfldy,bfldz)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param
      use corgan_com_M
      use cindex_com_M
      use numpar_com_M
      use cprplt_com_M
      use cophys_com_M
      use ctemp_com_M
      use blcom_com_M
!
      implicit none
!-----------------------------------------------
!   G l o b a l   P a r a m e t e r s
!-----------------------------------------------
     real(double) :: harris,flvx,flvy,flvz,bfldx,bfldy,bfldz,vt_shape
     integer :: np,kr
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------

!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
     integer :: mode
     real(double)::xc,zc,yc,wid,xcp,zcp,r,th,cth,sth, &
                   wspx,wspy,wspz,bth,bth0,jy,segno, pressure
!-----------------------------------------------
!   E x t e r n a l   F u n c t i o n s
!-----------------------------------------------
!
segno = sign(1.,qom(icoi(kr)))
!
!    Compute grid quantities
!
!
      ipjk = ijk + iwid
      ipjpk = ijk + iwid + jwid
      ijpk = ijk + jwid
!
      ijkp = ijk + kwid
      ipjkp = ijk + iwid + kwid
      ijpkp = ijk + jwid + kwid
      ipjpkp = ijk + iwid + jwid + kwid
!

      wspx = 0.125d0*(x(ipjk)+x(ipjpk)+x(ijpk)+x(ijk)+x(ipjkp)+x(ipjpkp)+x(ijpkp)+x(ijkp))
      wspy = 0.125d0*(y(ipjk)+y(ipjpk)+y(ijpk)+y(ijk)+y(ipjkp)+y(ipjpkp)+y(ijpkp)+y(ijkp))
      wspz = 0.125d0*(z(ipjk)+z(ipjpk)+z(ijpk)+z(ijk)+z(ipjkp)+z(ipjpkp)+z(ijpkp)+z(ijkp))

    xc=ibar*dx/2.d0
    yc=jbar*dy/2.d0
    zc=kbar*dz/2.d0

    wid=xc/10.d0*2
    mode=1
    xcp=xc+0d0*wid/2000.d0*sin(2.d0*pi*wspy*mode/dy/jbar)
    zcp=zc+0d0*wid/2000.d0*cos(2.d0*pi*wspy*mode/dy/jbar)

    r=sqrt((wspx-xcp)**2+(wspz-zcp)**2)
    th=acos((wspx-xcp)/(r+1.d-10))
    cth=(wspx-xcp)/(r+1.d-10)
    sth=(wspz-zcp)/(r+1.d-10)

    bth0=bx_ini
    bth=bth0*r/(r**2+wid**2+1.d-10)
    jy=bth0*2*wid**2/(r**2+wid**2+1.d-10)**2
    pressure=.5d0*bth0**2*wid**2/(r**2+wid**2+1.d-10)**2
    pressure=pressure+.001d0/4d0

    vt_shape=sqrt(pressure/(siepx(1)**2/abs(qom(1))+siepx(2)**2/abs(qom(2))))
    harris=1d0

!    rgeom=sqrt((wspx-xc)**2+(wspz-zc)**2)
!        if(rgeom.ge.5d0-min(dx,dz)/2d0) then
!        rhoprof=rhoprof*1d2
!        uprof=0d0
!        endif
! campo magnetico

    bfldy=by_ini
    bfldx=-sth*bth
    bfldz=cth*bth

    flvx=0d0
    if(segno.lt.0d0) then
    flvy=-segno*jy/harris
    else
    flvy=0d0
    end if
    flvz=0d0
!write(*,*) 'bbb',bfldx,bfldy,bfldz,flvy
end subroutine rope_equilibrium


subroutine tokamak_equilibrium(np,kr,harris,t_shape,flvx,flvy,flvz,bfldx,bfldy,bfldz)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param
      use corgan_com_M
      use cindex_com_M
      use numpar_com_M
      use cprplt_com_M
      use cophys_com_M
      use ctemp_com_M
      use blcom_com_M
!
      implicit none
!-----------------------------------------------
!   G l o b a l   P a r a m e t e r s
!-----------------------------------------------
     real(double) :: harris,flvx,flvy,flvz,bfldx,bfldy,bfldz,t_shape
     integer :: np,kr
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
!-----------------------------------------------
!   E x t e r n a l   F u n c t i o n s
!-----------------------------------------------
!
!   obsolete not used
!
!              particle part
!
!       Particle thermal velocity obtained from
!       siep (from input file) and T wall
!     siep ... core temperature
!     T_wall ... wall temperature
!     linear variation with zeta (k)
!
!       gradient to study ITG mode
!
!                     variation = (pzta(np)-2.)/float(kbar)
!                     scale=siep(kr)*(1.-variation)+T_wall(kr)*variation
!
!       Density Gradient from Birn 2D equilibrium for magnetic
!       reconnection in a planar current sheet.
!       empirical BIRN profile OR
!       x-dependence from a quasi parabolic form (Pritchett)
!

!   magnetic field part

!
!         xc=xc*toroid+strait
!         rmr = 1./sqrt(xc*xc+yc*yc)
!         sphi = yc*rmr
!         cphi = xc*rmr
!         stheta = ws/(wsr+1.e-35)
!         ctheta = zc/(wsr+1.e-35)
!
!     calculate toroidal magnetic field
!     Ref: Bauer, Betancourt, and Garabedian
!     "A Computational Method in Plasma Physics"
!
!
!      wsq0=q0*float(ibar)/float(jbar)
!      denom=2.*pi/(float(ibar)*float(kbar))
!
!      bxn(ijk)=btorus*(nuy(ijk)*(wsq0*etaz(ijk)+tsiz(ijk))
!     &             -nuz(ijk)*(wsq0*etay(ijk)+tsiy(ijk)))*denom
!     &     *(rmaj/wsmr+strait)*wsr
!
!      byn(ijk)=btorus*(nuz(ijk)*(wsq0*etax(ijk)+tsix(ijk))
!     &             -nux(ijk)*(wsq0*etaz(ijk)+tsiz(ijk)))*denom
!     &     *(rmaj/wsmr+strait)*wsr
!
!      bzn(ijk)=btorus*(nux(ijk)*(wsq0*etay(ijk)+tsiy(ijk))
!     &             -nuy(ijk)*(wsq0*etax(ijk)+tsix(ijk)))*denom
!     &     *(rmaj/wsmr+strait)*wsr
!     &     +bvertcl
!
!      wsnr=(k-1.5)/(kbp1-1.)
!      phipot(ijk)=0.5*(1.-wsnr**2)*efld0
!

end subroutine tokamak_equilibrium

subroutine boundary_state(ib,jb,kb,density,temperature,flvx,flvy,flvz,bfldx,bfldy,bfldz)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      use vast_kind_param, ONLY:  double
      use corgan_com_M
      use cindex_com_M
      use numpar_com_M
      use cprplt_com_M
      use cophys_com_M
      use ctemp_com_M
      use blcom_com_M
!
      implicit none
!-----------------------------------------------
!   G l o b a l   P a r a m e t e r s
!-----------------------------------------------
    integer :: ib,jb,kb
    real(double) :: density,temperature,flvx,flvy,flvz,bfldx,bfldy,bfldz
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
    integer :: np, kr
    real(double) :: harris, t_shape,efldx,efldy,efldz
!-----------------------------------------------
!   E x t e r n a l   F u n c t i o n s
!!-----------------------------------------------
!-----------------------------------------------

    ijk= 1 + (ib - 1)*iwid + (jb - 1)*jwid + (kb - 1)*kwid
    np=1

    density=0d0
    temperature=0d0
    do kr=1,nrg

    call initial_state(np,kr,harris,t_shape,flvx,flvy,flvz,bfldx,bfldy,bfldz,efldx,efldy,efldz)
    density=density+harris*rhr(kr)
    temperature=temperature+t_shape

    enddo

    temperature=temperature*.25d0
    density=density*.5d0

end subroutine boundary_state

subroutine boundary_state_ijk(density,temperature,flvx,flvy,flvz,bfldx,bfldy,bfldz)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      use vast_kind_param, ONLY:  double
      use corgan_com_M
      use cindex_com_M
      use numpar_com_M
      use cprplt_com_M
      use cophys_com_M
      use ctemp_com_M
      use blcom_com_M
!
      implicit none
!-----------------------------------------------
!   G l o b a l   P a r a m e t e r s
!-----------------------------------------------
    real(double) :: density,temperature,flvx,flvy,flvz,bfldx,bfldy,bfldz
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
    integer :: np, kr
    real(double) :: harris, t_shape,efldx,efldy,efldz
!-----------------------------------------------
!   E x t e r n a l   F u n c t i o n s
!!-----------------------------------------------
!-----------------------------------------------

    np=1

    density=0d0
    temperature=0d0
    do kr=1,nrg

    call initial_state(np,kr,harris,t_shape,flvx,flvy,flvz,bfldx,bfldy,bfldz,efldx,efldy,efldz)
    density=density+harris*rhr(kr)
    temperature=temperature+t_shape

    enddo

    temperature=temperature/dble(nrg)
    density=density/dble(nsp)

end subroutine boundary_state_ijk

subroutine onewave_equilibrium(np,kr,den_shape,t_shape,flvx,flvy,flvz,bfldx,bfldy,bfldz,efldx,efldy,efldz)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param
      use corgan_com_M
      use cindex_com_M
      use numpar_com_M
      use cprplt_com_M
      use cophys_com_M
      use ctemp_com_M
      use blcom_com_M
!
      implicit none
!-----------------------------------------------
!   G l o b a l   P a r a m e t e r s
!-----------------------------------------------
     real(double) :: den_shape,flvx,flvy,flvz,bfldx,bfldy,bfldz,t_shape,efldx,efldy,efldz
     integer :: np,kr
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
	real(double) :: xc, zc
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
!-----------------------------------------------
!   E x t e r n a l   F u n c t i o n s
!-----------------------------------------------
!-----------------------------------------------!
!    Compute grid quantities
!

                     ipjk = ijk + iwid
                     ipjpk = ijk + iwid + jwid
                     ijpk = ijk + jwid
!
                     ijkp = ijk + kwid
                     ipjkp = ijk + iwid + kwid
                     ijpkp = ijk + jwid + kwid
                     ipjpkp = ijk + iwid + jwid + kwid
!
                     xc = 0.125d0*(x(ipjk)+x(ipjpk)+x(ijpk)+x(ijk)+x(ipjkp)+x(ipjpkp)+x&
                          (ijpkp)+x(ijkp))
                     zc = 0.125d0*(z(ipjk)+z(ipjpk)+z(ijpk)+z(ijk)+z(ipjkp)+z(ipjpkp)+z&
                          (ijpkp)+z(ijkp))
!! Unifrom density
den_shape=1d0
! Uniform temperature
t_shape=1.d0
!
! Plasma flow
!


if(kbar>ibar) then
flvy=vvi(kr)- wave_amplitude *cos(zc/kbar/dz*2d0*pi)
flvx = uvi(kr) !- wave_amplitude *cos(zc/kbar/dz*2d0*pi)
flvz = wvi(kr)
else
flvy=vvi(kr)- wave_amplitude *sin(xc/ibar/dx*2d0*pi)
flvz = wvi(kr) ! - wave_amplitude *sin(xc/ibar/dx*2d0*pi)
flvx = uvi(kr)
end if
!
! Magnetic Field
!     bfldx=bx_ini*tanh((zc-.5*(zf+ze))/el_rec)
bfldx = bx_ini
if(kbar>ibar) then
bfldy = by_ini +wave_amplitude *cos(zc/kbar/dz*2d0*pi)
else
bfldy = by_ini +wave_amplitude *sin(xc/ibar/dx*2d0*pi)
end if
bfldz = bz_ini

 efldx = -( flvy*bfldz-flvz*bfldy)
 efldy = -(-flvx*bfldz+flvz*bfldx)
 efldz = -( flvx*bfldy+flvy*bfldx)

end subroutine onewave_equilibrium


subroutine doubletearing_equilibrium(np,kr,harris,t_shape,flvx,flvy,flvz,bfldx,bfldy,bfldz)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param
      use corgan_com_M
      use cindex_com_M
      use numpar_com_M
      use cprplt_com_M
      use cophys_com_M
      use ctemp_com_M
      use blcom_com_M
!
      implicit none
!-----------------------------------------------
!   G l o b a l   P a r a m e t e r s
!-----------------------------------------------
     real(double) :: harris,flvx,flvy,flvz,bfldx,bfldy,bfldz,t_shape
     integer :: np,kr
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
   real(double) :: fx,fxone,fxtwo,shaperx,shapery,segno,xc,zc,xone,xtwo,zone,ztwo,zzc
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
!-----------------------------------------------
!   E x t e r n a l   F u n c t i o n s
!-----------------------------------------------
!-----------------------------------------------!
!    Compute grid quantities
!
!
                     ipjk = ijk + iwid
                     ipjpk = ijk + iwid + jwid
                     ijpk = ijk + jwid
!
                     ijkp = ijk + kwid
                     ipjkp = ijk + iwid + kwid
                     ijpkp = ijk + jwid + kwid
                     ipjpkp = ijk + iwid + jwid + kwid
!

                     xc = 0.125d0*(x(ipjk)+x(ipjpk)+x(ijpk)+x(ijk)+x(ipjkp)+x(ipjpkp)+x&
                          (ijpkp)+x(ijkp))
                     zc = 0.125d0*(z(ipjk)+z(ipjpk)+z(ijpk)+z(ijk)+z(ipjkp)+z(ipjpkp)+z&
                          (ijpkp)+z(ijkp))
!
                     zone = 0.25d0*(z(ipjk)+z(ipjpk)+z(ijpk)+z(ijk))
                     ztwo = 0.25d0*(z(ipjkp)+z(ipjpkp)+z(ijpkp)+z(ijkp))
                     xtwo = 0.25d0*(x(ipjk)+x(ipjpk)+x(ipjkp)+x(ipjpkp))
                     xone = 0.25d0*(x(ijk)+x(ijpk)+x(ijkp)+x(ijpkp))
!
   harris = 1d0!.1d0+1.d0/cosh((zc+kbar*dz/4d0)/el_rec)**2/el_rec+1.d0/cosh((zc-kbar*dz/4d0)/el_rec)**2/el_rec
t_shape=1d0
fx=1.d0
fxone=1d0
fxtwo=1d0
segno = sign(1.,qom(icoi(kr)))
shaperx = tanh((pz(np)-0.5*(zf+ze))/el_rec)/cosh((pz(np)-0.5*(zf+ze))/el_rec)/el_rec
shapery = 1.d0/cosh((pz(np)-0.5*(zf+ze))/el_rec)**2/el_rec
!
flvx=0
flvy=0
flvz=0
! uvi,vvi must be chosen to support the current
! campo magnetico

!bfldx = bx_ini*sin(2d0*pi*((zc-kbar*dz/4d0)/kbar/dz))
zzc=zc-kbar*dz/2d0
bfldx = bx_ini*(tanh((zzc+kbar*dz/4d0)/el_rec)-tanh((zzc-kbar*dz/4d0)/el_rec)-1d0)
bfldy = by_ini/cosh((zzc+kbar*dz/4d0)/el_rec)-by_ini/cosh((zzc-kbar*dz/4d0)/el_rec)
bfldz = bz_ini

end subroutine doubletearing_equilibrium
